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Abstract 

A  comprehensive,  steady-state,  computational  model  of  a  proton  exchange  membrane  fuel  cell  (PEMFC)  derived  from  first  principles  is 
presented.  The  model  is  two-dimensional  and  includes  the  transport  of  liquid  water  within  the  porous  electrodes  as  well  as  the  transport  of 
gaseous  species,  protons,  energy,  and  water  dissolved  in  the  ion  conducting  polymer.  Electrochemical  kinetics  are  modeled  with  standard 
rate  equations  adapted  to  an  agglomerate  catalyst  layer  structure.  Some  of  the  physical  properties  used  in  constructing  the  model  are 
determined  experimentally  for  an  in-house  membrane  electrode  assembly  (ME A)  and  are  presented  herein.  Experimental  results  obtained 
for  the  MEA  are  used  to  validate  the  computational  model.  Modeling  results  are  presented  that  illustrate  the  importance  of  the  transport  of 
water  within  the  porous  sections  of  the  cell  and  in  the  polymer  regions  of  the  MEA. 
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1.  Introduction 

The  development  of  economical,  commercially  available 
proton  exchange  membrane  (PEM)  fuel  cells  is  essential  to 
the  development  of  a  hydrogen  based  energy  infrastructure. 
The  PEM  fuel  cell  is  particularly  important  to  the  develop¬ 
ment  of  hydrogen  fueled  vehicles  as  demonstrated  by  the  fact 
that  all  of  the  major  automobile  manufacturers  are  currently 
involved  in  some  level  of  PEM  fuel  cell  research  with  the 
goal  of  widespread  commercialization  within  10-15  years. 
To  realize  this  and  other  development  goals,  fuel  cell  tech¬ 
nology  must  continue  to  improve  with  respect  to  cost  and 
performance.  Computational  models  of  fuel  cells  can  con¬ 
tribute  to  this  effort  by  providing  researchers  with  the  ability 
to  create  and  optimize  fuel  cell  designs  rapidly  and  inexpen¬ 
sively  before  actually  building  a  prototype. 

1.1.  Overview  of  fuel  cell  models 

The  need  for  computational  tools  to  support  the  design 
process  has  led  to  the  development  of  a  number  of  fuel  cell 
models.  These  models  can  generally  be  characterized  by  the 
scope  of  the  model.  In  many  cases,  modeling  efforts  focus  on 
a  specific  part  or  parts  of  the  fuel  cell,  like  the  cathode  cat¬ 
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alyst  layer  [1-4],  the  cathode  electrode  (gas  diffusion  layer 
plus  catalyst  layer)  [5],  or  the  membrane  electrode  assem¬ 
bly  (MEA)  [6,7].  These  models  are  very  useful  in  that  they 
may  include  a  large  portion  of  the  relevant  fuel  cell  physics 
while  at  the  same  time  having  relatively  short  solution  times. 
However,  these  narrowly  focused  models  neglect  important 
parts  of  the  fuel  cell  making  it  impossible  to  get  a  com¬ 
plete  picture  of  the  phenomena  governing  fuel  cell  behav¬ 
ior.  Models  that  include  all  parts  of  a  fuel  cell  are  typically 
two-  or  three-dimensional  and  reflect  many  of  the  physical 
processes  occurring  within  the  fuel  cell  [8-12].  Comprehen¬ 
sive  models  rely  on  the  determination  of  a  large  number  of 
properties  and  operating  parameters  and  can  be  much  more 
computationally  intensive,  leading  to  longer  solution  times. 
However,  these  disadvantages  are  typically  outweighed  by 
the  benefit  of  being  able  to  assess  the  influence  of  a  greater 
number  of  design  parameters  and  their  associated  physical 
processes.  The  model  presented  in  this  work  is  a  compre¬ 
hensive  two-dimensional  model  that  includes  multicompo¬ 
nent  and  multiphase  transport  both  along  the  flow  direction 
(down  the  gas  channel)  and  through  the  MEA. 

1.2.  Models  of  liquid  water  transport 

A  comprehensive  computational  model  should  include 
the  equations  and  other  numerical  relations  needed  to  fully 
define  fuel  cell  behavior  over  the  range  of  interest.  Early 
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Nomenclature 

Parameters  and  variables 
c  concentration  (mol/mm3) 

ceff  effective  heat  capacity  (J/g  K) 

D  diffusivion  coefficient  (mm2/s) 

E{ h  Nernst  potential  (V) 

F  Farraday’s  constant 

hfg  enthalpy  of  vaporization  (J/g) 

&eff  effective  thermal  conductivity  (W/mm  K) 

M  molar  mass  (g/mol) 

n  number  of  electrons 

P  pressure  (Pa) 

R  universal  gas  constant  (J/mol  K) 

.v  saturation 

standard  entropy  of  formation  (J/mol  K) 

S  source  term — see  Table  2 

t  thickness  (mm) 

T  temperature  (K) 

u  velocity  (mm/s) 

Veen  cell  voltage  (V) 

w  mass  fraction 

Greek  letters 

<fi  potential  (V) 

£  volume  fraction 

y  switch  function 

k  relative  humidity 

A  polymer  water  content  H2  O/S 03“ 

pc  viscosity  (Pas) 

p  density  (g/mm3) 

a  ionic  conductivity  (fimm)-1 

r  tortuosity 

Subscripts 

a  anode 

act  activation 

c  cathode 

C  carbon 

coll  current  collector 

d  index  for  electrodes  (anode:  d  =  a, 

cathode:  d  =  c) 
drag  electro-osmotic  drag 

Dar  Darcy  pressure  loss 

eff  effective 

evap  evaporation 

H2  hydrogen 

i  ionic 

LV  mass  transfer  from  liquid  to  vapor 

O2  oxygen 

p  polymer 

pi  plate 

rev  reversible  heat 

v  vapor 

void  void  space 


W  water 

WD  water  dissolved  in  polymer 

WV  water  vapor 

WL  water  liquid 

WP  water  production 

v  v-direction  -  through  MEA 

y  y-direction  -  down  channel 

Q  Ohmic 

Superscripts 

cat  catalyst  layer 

cp  capillary  pores 

g  gas  phase 

gdl  gas  diffusion  layer 

p  polymer  phase 

v  vapor  phase 


multidimensional  models  described  gas  transport  in  the  flow 
channels,  gas  diffusion  layers  (GDLs),  and  the  membrane 
[11,10].  More  recent  models  include  a  detailed  description 
of  the  catalyst  layers  that  reflects  the  agglomerate  nature  of 
these  regions  [12].  Recently,  there  has  been  an  interest  in  de¬ 
scribing  operating  regimes  that  are  dominated  by  mass  trans¬ 
port  limitations  resulting,  in  part,  from  the  formation  and 
transport  of  liquid  water  within  the  fuel  cell.  To  model  fuel 
cell  performance  in  these  regimes,  it  is  necessary  to  include 
equations  that  describe  not  only  the  motion  of  water  within 
the  liquid  phase,  but  also  mass  transfer  between  phases.  He 
et  al.  [15]  present  a  model  for  liquid  water  transport  in  the 
cathode  gas  diffusion  layer  of  a  fuel  cell  with  an  interdigi- 
tated  flowfield.  As  part  of  their  work  they  formulate  a  switch¬ 
ing  function  that  is  used  to  toggle  source  terms  on  and  off  de¬ 
pending  on  whether  water  is  condensing  into  the  liquid  phase 
or  evaporating  into  the  vapor  phase.  They  also  introduce  a 
term  to  account  for  water  transport  by  advection,  or  transport 
due  to  the  motion  of  the  bulk  flow,  which  is  an  important 
transport  mechanism  for  interdigitated  flowfields.  In  a  con¬ 
tinuation  of  this  work,  Natarajan  and  Nguyen  [16]  develop  a 
diffusive  expression  to  account  for  water  transport  via  cap¬ 
illary  pressure  in  the  porous  electrode.  They  neglect  water 
transport  by  advection  since  their  model  uses  a  conventional 
flowfield.  Wang  and  Wang  [17]  have  also  developed  a  model 
for  two  phase  flow  in  an  interdigitated  cathode.  They  as¬ 
sume  transport  by  capillary  pressure  only  and  do  not  include 
mass  transport  between  phases.  All  three  of  these  models 
treat  the  catalyst  layer  as  an  interface  and  so  do  not  con¬ 
sider  effects  due  to  water  buildup  within  the  catalyst  layer. 
The  model  of  Baschuk  and  Li  [18]  is  one-dimensional  and 
includes  liquid  water  effects  in  both  the  gas  diffusion  layer 
and  the  catalyst  layer.  In  the  Baschuk  and  Li  model,  phase 
change  is  neglected  and  the  volume  fraction  of  liquid  water 
in  the  porous  regions  must  be  specified  as  opposed  to  being 
calculated. 
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1.3.  Determination  of  physical  properties 

Another  important  part  of  developing  a  useful  computa¬ 
tional  fuel  cell  model  is  the  accurate  determination  of  the 
physical  properties  on  which  the  model  is  based.  These 
properties  may  include:  porosity  of  the  gas  diffusion  layer 
and  catalyst  layer;  the  ionic  conductivity  of  the  polymer 
membrane;  mass  transfer  coefficients;  kinetic  parameters 
such  as  the  exchange  current  density  and  reaction  surface 
area;  and  structural  properties  such  as  the  thicknesses  of  the 
catalyst  layers,  membrane,  and  gas  diffusion  layers.  Many 
of  these  properties  can  be  experimentally  determined  with 
considerable  accuracy.  Others,  like  interphase  mass  transfer 
coefficients,  are  more  difficult  to  determine  and  must  be 
estimated.  In  either  case,  it  is  important  to  take  great  care 
with  the  values  used  as  they  may  significantly  affect  the 
results  [1]. 

Prior  research  has  determined  many  of  the  physical  prop¬ 
erties  needed  for  modeling.  The  work  by  Parthasarathy  et  al. 
[13]  and  Zhang  et  al.  [19]  contains  experimentally  deter¬ 
mined  values  for  many  of  the  kinetic  parameters  involved 
with  the  oxygen  reduction  reaction  occurring  at  the  cathode 
as  well  as  data  regarding  the  solubility  and  diffusivity  of 
oxygen  in  Nation®.  Ihonen  et  al.  [14]  present  data  for  the 
pore  size  distribution  and  porosity  of  the  catalyst  layer.  In 
their  work,  they  compare  the  pore  size  distribution  as  de¬ 
termined  by  both  gas  and  mercury  porosimetry  and  show 
that  the  two  methods  provide  nearly  identical  results  for  the 
catalyst  layers  that  they  tested.  Many  of  the  required  prop¬ 
erties  for  Nation®  1100,  including  ionic  conductivity  and 
the  electro-osmotic  drag  coefficient,  can  be  found  in  the 
work  by  Springer  et  al.  [6].  Marr  and  Li  [2]  present  data 
for  catalyst  surface  area  as  it  relates  to  the  type  of  catalyst 
used  in  the  ME  A.  This  data  is  for  a  perfectly  uniform  cata¬ 
lyst  deposition.  However,  transmission  electron  microscope 
(TEM)  images  such  as  those  shown  by  Siegel  et  al.  [12]  re¬ 
veal  that  the  catalyst  deposition  is  not  uniform  and  property 
data  must  be  modified  to  account  for  the  effect  of  catalyst 
non-uniformities.  In  the  current  work,  we  experimentally 
determine  the  catalyst  layer  porosity  and  thickness  using  a 
scanning  electron  microscope  (SEM)  and  the  reaction  sur¬ 
face  area  using  cyclic  voltammetry. 

1.4.  Characteristics  of  present  model 

The  present  work  presents  a  comprehensive  model  of 
a  PEM  fuel  cell  that  incorporates  the  significant  physical 
processes  and  the  key  parameters  affecting  fuel  cell  per¬ 
formance.  The  model  is  two-dimensional  and  includes  the 
transport  of  gaseous  species,  protons,  energy,  and  water 
dissolved  in  the  ion  conducting  polymer.  The  model  also 
addresses  the  transport  of  liquid  water  within  the  fuel  cell 
with  liquid  water  assumed  to  be  transported  by  capillary 
pressure  within  the  gas  diffusion  layers  and  catalyst  layers 
and  by  advection  within  the  gas  channels.  Water  is  assumed 
to  be  exchanged  among  three  phases — liquid,  vapor,  and 


dissolved1 — and  equilibrium  among  these  phases  is  as¬ 
sumed.  Electrochemical  kinetics  are  modeled  with  standard 
rate  equations  adapted  to  an  agglomerate  catalyst  structure. 
The  model  reflects  the  influence  of  numerous  parameters 
on  fuel  cell  performance  including  geometry,  porosity  of 
the  cell  materials,  catalyst  area,  polymer  properties,  cata¬ 
lyst  layer  composition,  and  others.  This  paper  describes  the 
development  of  the  model,  the  determination  of  properties 
for  use  in  the  model,  the  validation  of  the  model  using  ex¬ 
perimental  data,  and  the  application  of  the  model  to  explain 
observed  experimental  phenomena. 

2.  Model  development 

Fig.  1  shows  the  solution  domain  of  the  model.  In  the 
anode  and  cathode  gas  channels,  fuel  and  oxidant  flow 
along  the  surface  of  the  membrane  electrode  assembly.  In 
these  regions,  the  flow  is  considered  to  be  laminar.  Reac¬ 
tants  move  from  the  gas  channels  into  the  gas  diffusion 
layers  (GDL)  which  serve  to  more  uniformly  distribute  the 
reactants  across  the  surface  of  the  catalyst  layer.  In  the  cat¬ 
alyst  layers,  the  reactants  are  transported  by  diffusion  and 
advection  while  participating  in  electrochemical  reactions. 
The  polymer  membrane  is  assumed  to  transport  only  pro¬ 
tons  and  dissolved  water.  The  current  collector  plates  are 
not  explicitly  included  in  the  computational  model.  Their 
influence  with  respect  to  energy  transport  is  included  and 
discussed  in  the  next  section. 

2.1.  Governing  equations 

The  governing  equations  include  conservation  of  mass, 
momentum,  ionic  charge,  and  energy  as  well  as  individual 
species  conservation  equations.  The  governing  equations  are 
presented  in  Table  1.  Related  terms  that  describe  the  rate 
at  which  the  conserved  quantity  is  added  or  removed  from 
the  solution  domain  are  referred  to  as  “source”  terms  and 
are  presented  in  Table  2.  Constitutive  relations  describing 
reaction  rates,  polymer  properties,  gas  properties,  and  liquid 
properties  are  provided  in  Table  3. 

Each  of  the  governing  equations  is  developed  for  the 
entire  solution  domain  of  the  model  with  source  and  trans¬ 
port  terms  modified  in  each  region  to  reflect  the  appropriate 
physical  phenomena.  Eq.  (1)  describes  conservation  of  mass 
for  the  entire  gas  phase.  Source  terms  reflect  changes  in  the 
overall  gas  phase  mass  flow  due  to  consumption  or  produc¬ 
tion  of  gas  species  via  reaction  and  due  to  mass  transfer 
between  water  in  the  vapor  phase  and  water  that  is  in  the  liq¬ 
uid  phase  or  dissolved  in  the  polymer.  The  system  of  equa¬ 
tions  that  model  the  interphase  mass  transfer  is  described 
in  more  detail  in  the  following  section.  The  gas  phase  mix¬ 
ture  density  is  determined  from  the  temperature,  pressure, 


1  Here,  “dissolved”  refers  to  water  that  exists  within  the  polymer  mem¬ 
brane  or  the  polymer  phase  of  the  catalyst  layer. 
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Table  1 

Governing  equations 
Conservation  equation 
Mass 

Momentum 

Oxygen 
Water  vapor 

Liquid  water 

Dissolved  water 

Ionic  charge 
Thermal  energy 


Fig.  1.  Solution  domain. 


General  form  of  equation 


V  •  ( pgu )  =  Sho  +  Sq2  +  Slv  +  SwpKlv  —  SwdKwd 


U  •  —  T  pF7  Ux  T  6f)ar5x 

dx 

dP  2 

U  •  V[p  Uy\  —  -  +  IxW  Uy  +  ^Dar,y 

dy 

pgu  ■  Wwo2  +  wo2W  ■  (pgu)  =  V  •  (£>52pgV wo2)  +  So2 

pgu  ■  Vinwv  +  rPwvV  •  (pgu)  =  V  •  (D^vpgVu;wv)  +  Slv  +  SwpKlv  —  SwdKwd 


_  cp  ^  .  Slv  ,  ^wp(1  —  Klv)  6wd0  -  Kwd) 

u  -  Vs  =  V  •  (Df)LWs) - I - 

PWL  PWL  PWL 

v  •  (XepOiV*)  =  V  ■  (Z4D  v4d)  +  ^ 

V  •  (£p(JiV0i)  +  S[  =  0 


V  •  (XspC'i V0i)  =  V  •  (D^d  VcP,d)  + 


pgceffM  •  vr  =  v  •  (Peff  vr)  +  +  6*rcv  T  5aCt  T  Spc  +  Sp\ 


(1) 

(2a) 

(2b) 

(3) 

(4) 

(5) 

(6) 

(V) 

(8) 


Table  2 
Source  terms 

Source  term  (region  of  application  —  zero  in  other  regions) 

Darcy  pressure  drop  in  the  jc-direction  (Pa/mm)  (anode  and  cathode  GDLs  and  catalyst  layers) 
Darcy  pressure  drop  in  the  y-direction  (Pa/mm)  (anode  and  cathode  GDLs  and  catalyst  layers) 

Hydrogen  reaction  rate  (g/mm3  s)  (anode  catalyst  layer) 

Oxygen  reaction  rate  (g/mm3  s)  (cathode  catalyst  layer) 

Water  production  rate  (g/mm3  s)  (cathode  catalyst  layer) 

Mass  transfer  rate  from  liquid  to  vapor  (g/mm3  s)  (anode  and  cathode  cat.  layers,  GDLs, 
gas  channels) 

Mass  transfer  rate  into  the  dissolved  phase  (g/mm3  s)  (anode  and  cathode  catalyst  layers) 

Ionic  current  generation  (A/mm3)  (anode  and  cathode  catalyst  layers) 

Heat  source  due  to  ohmic  heating,  (anode  and  cathode  catalyst  layers;  membrane) 

Heat  source  due  to  reversible  chemical  reaction  (anode  and  cathode  catalyst  layers) 

Heat  source  due  to  activation  loss  (W/mm3)  (anode  and  cathode  catalyst  layers) 

Heat  source  due  to  phase  change  (W/mm3)  (cathode  catalyst  layer) 


Defining  equation 


^Dar,*  —  ux 

K 

^Dar,y  —  Uy 

K  " 


SHi  = 


SCh  = 


MH2 

IF 

Mo2 


|Reff,< 


|Reff,< 


Mw 

Swp  -  -^rlReff.cl 

Sl\  =  fsyix  -  f(\  -  s)(l  -  Klv) 


3\vD  —  bm  (PwV  Pwv) 


5)  —  Reff.d 


+  at  anode  (d  =  a) 
—  at  anode  (d  —  c) 


Sq  =  V0j  •  (ep(jiV0i) 
o  Poff,kT\~^n 


Srpv  — 


nkF 


5act  —  (0e,d  0i)Reff,A' 

6pc  =  (~Slv  +  3’dw  Kdw  )b{'g 

g  _  (TColl  —  T)  /  ^gdl  ^  tcoll 

^gdl  \  ^gdl  ^coll 


(9) 

(10) 

(ID 

(12) 

(13) 

(14) 

(15) 

(16) 

(17) 

(18) 

(19) 

(20) 


Heat  source  due  to  the  current  collector  plate  (W/mm3)  (anode  and  cathode  GDLs) 


-l 


(21) 


Table  3 

Closure  relations 


N.P  Siegel  et  al.  /  Journal  of  Power  Sources  128  (2004)  173-184 


111 


Parameter 


Defining  equation 


Reaction  rate  equations  (anode:  d  =  a,  k  =  H2,  = 

Effective  reaction  rate,  Reff  (A/mm3) 

Reaction  rate,  R  (A/mm3) 


1;  cathode:  d  =  c,  k  =  O2,  no,  —  1) 

^eff,d  —  /^d^aggl,d 

R&  —  (1  V)Apvz0,k 


Y 


d'k.rcf 


^e(0e,d-0i)o'd(«di7/^7)  _  e-(0e,d-0i)°'d(«di7/^:OJ 


Electrical  potential,  0C 
Agglomerate  effectiveness,  r/agg\ 


0e,a  —  0 

0e,c  —  (^th  Well) 
_  3  /  1 

"aggl’d  “  A,  \tanh(^d) 


1 

ft 


Thiele’s  modulus,  f> 

Dissolved  gas  concentration,  cp  (mol/mm3) 

Material  properties 

Polymer  water  content,  k  (mol-W/mol-p) 

Ionic  conductivity  of  the  polymer,  cri  (S/mm) 

Water  diffusivity  in  the  polymer,  D^w  (mm2/s) 

Water  vapor  activity,  a 

Density  of  water  vapor  in  equilibrium  with  the  polymer, 
Pwvp  (g/mm3) 

Gas  density,  pg  (g/mm3) 


o  j  ,  l^dl 
Pd  ~  ^ J  c{D{ndF 


cp  _  hPcsT 

ck  —  nkckJ 


k 


CPWMP 

Pm 


ori  =  (0.0005139  k  -  0.000326)  exp 

1 


D^w  =  1.3  x  10  4  exp 


a  =  1.76e-6A.4  +  2.17  e-4/G  —  8.80e_JA.z  +  0.16  A.  —  0.12 


2416 


-4  1  3 


1268.0 


1 


1 


1 


303  T 


303  T 


3i2 


%V  ~  ASAT  * 


a 


Pg  - 


RT  E (ujj/Mj) 


Gas  component  diffusivity  (mm2/s) 


Capillary  diffusivity,  D^L  (mm2/s) 


Effective  thermal  conductivity  of  region  r  (W/mmK) 


D 

D 


g  =  Pg(  1  -yk) 
E/^/c  (y.j/  Djk) 

Dgsgdl  a  - 
g  _  ^kWoicfl1 

eff,k  fgdl 


£)§  _  /)grpgdl  n  — 

^eff.k  — 


kL  void 


gas  channel 
5) 

—  porous  GDL 
s)] 1  5  porous  catalyst  layer 


PWLg 

Mwl 


K(s) 


dlf 

ds 


dP 

K(s) — -  =  0.0155/  -  0.0213/  +  0.0088s  +  0.0002 
ds 


Wff  =  CidC1  -  s)kg  +  Cid^WL  +  £p^P  +  £c^c 


Equilibrium  control  functions 

Pg 

Liquid/vapor  switch,  klv  Klv  — ►  1  for  <  0.98,  elseO 

Psat 

Dissolved  water  switch,  kwd  Kwd  =  lforp^y  <  p|yV,  elseO 


Klv  =  1—0.5 


1  +  tanh  |  61  -59 

PSAT 


Kwd  =  0.5  + 


%v  %v 


Pwv  Pwv 


(22) 

(23) 

(24) 

(25) 

(26) 

(27) 

(28) 

(29) 

(30) 

(31) 

(32) 

(33) 

(34a) 

(34b) 

(34c) 

(35a) 

(35b) 

(36) 

(37) 

(38) 


and  mixture  molar  mass  using  the  ideal  gas  equation, 
Eq.  (33). 

Conservation  of  momentum  is  expressed  by  the 
Navier-Stokes  equations  in  vector  form,  Eqs.  (2a)  and  (2b), 
which  are  modified  by  source  terms  described  by  Eqs.  (9) 
and  (10)  to  account  for  Darcy  flow  in  the  porous  regions 
of  the  model.  The  Darcy  terms  are  active  in  the  GDLs  and 
catalyst  layers  only;  the  inertial  and  bulk  viscous  terms  are 
neglected  in  these  regions. 

The  species  equations  for  oxygen  and  water  vapor  are 
given  by  Eqs.  (3)  and  (4),  respectively.  For  both  the  an¬ 
ode  and  cathode  sides  of  the  cell,  the  number  of  species 


equations  is  one  less  than  the  number  of  species.  At  the 
anode,  the  species  equation  for  water  vapor  is  solved.  The 
hydrogen  mass  fraction  is  calculated  from  the  solution  to 
the  water  vapor  equation  and  the  overall  gas  phase  con¬ 
servation  equation.  At  the  cathode,  both  water  vapor  and 
oxygen  species  equations  are  solved.  The  nitrogen  mass 
fraction  is  determined  from  the  solutions  to  these  species 
equations  and  the  overall  gas  phase  conservation  equa¬ 
tion.  The  diffusion  of  each  species  within  the  bulk  flow  is 
given  by  the  first  term  on  the  right  side.  The  diffusion  co¬ 
efficients  for  multi-component  flow,  which  are  determined 
by  Eqs.  (34a)-(34c),  are  based  on  a  simplification  of  the 


178 


N.P  Siegel  et  al.  /  Journal  of  Power  Sources  128  (2004)  173-184 


Stefan-Maxwell  equations  and  in  porous  regions  are  modi¬ 
fied  by  porosity  and  tortuosity  factors  [9,22].  The  advective 
term  for  each  of  the  species  equations  is  separated  into  two 
parts.  The  first  is  equal  to  the  product  of  velocity,  mixture 
density,  and  mass  fraction  gradient.  The  second  term  is  the 
product  of  mass  fraction  and  the  divergence  of  total  mass 
flux.  For  the  oxygen  species  equation,  Eq.  (3),  the  con¬ 
sumption  of  oxygen  via  reaction  is  reflected  in  the  source 
term  given  by  Eq.  (12).  The  water  vapor  transport  equation, 
Eq.  (4),  reflects  the  movement  of  water  vapor  by  diffusion 
and  bulk  motion;  the  exchange  of  water  between  the  vapor 
phase  and  the  liquid  phase  is  given  by  Eq.  (14);  the  ex¬ 
change  of  water  between  the  vapor  phase  and  the  dissolved 
phase  by  Eq.  (15);  and  the  production  of  water  by  Eq.  (13). 

Liquid  water  transport  within  the  cell  is  modeled  with 
Eq.  (5).  The  variable  that  describes  liquid  water  is  the  satura¬ 
tion,  s,  which  is  the  volume  fraction  of  liquid  water  relative 
to  the  pore  volume  in  the  porous  sections  of  the  fuel  cell,  i.e. 


fpore 

Within  the  porous  electrodes,  liquid  water  is  transported  by 
capillary  pressure  and  interphase  mass  transfer.  The  bulk 
flow  term  is  neglected.  The  diffusion  coefficient  for  liquid 
water,  D^L,  accounts  for  water  motion  via  capillary  pres¬ 
sure  and  is  based  on  a  semi-empirical  relation  between  cap¬ 
illary  pressure  and  saturation  [16].  Within  the  gas  channels, 
capillary  effects  are  neglected  and  the  liquid  is  assumed  to 
travel  as  droplets  of  negligible  volume  with  a  velocity  that  is 
equal  to  the  bulk  gas  velocity.  Water  is  exchanged  with  the 
liquid  phase,  due  to  evaporation  into  the  vapor,  as  given  by 
Eq.  (14);  production  from  the  cathode  reaction  as  given  by 
Eq.  (13)  if  the  adjacent  vapor  is  saturated;  and  transfer  into 
the  liquid  from  the  dissolved  phase  as  given  by  Eq.  (15). 
Equilibrium  between  the  various  water  phases  is  described 
in  more  detail  in  Section  2.2. 

Water  in  the  dissolved  form  is  transported  within  the  poly¬ 
mer  membrane  and  the  polymer  phase  of  the  catalyst  layer 
as  described  by  Eq.  (6).  The  transport  mechanisms  include 
diffusion,  electro-osmotic  drag,  and  inter-phase  mass  trans¬ 
fer  as  defined  by  Eq.  (15).  An  advective  term  would  be  re¬ 
quired  if  the  gas  pressure  in  the  anode  differed  significantly 
from  that  in  the  cathode.  For  this  model,  the  anode  and  cath¬ 
ode  pressures  are  assumed  to  be  equal  and  so  the  convec¬ 
tive  transport  of  water  through  the  MEA  is  neglected.  The 
transport  of  ions  in  the  polymer  regions  of  the  fuel  cell  is 
described  by  Eq.  (7),  where  the  source  term,  Eq.  (16),  rep¬ 
resents  the  production/consumption  of  protons  via  the  elec¬ 
trochemical  reactions  in  the  catalyst  layers. 

The  rate  of  the  electrochemical  reaction  is  described 
by  Eqs.  (22)  and  (23)  which  represent  a  form  of  the 
Butler- Volmer  relation  modified  by  an  effectiveness  factor 
to  account  for  the  effect  of  diffusion  through  the  agglom¬ 
erate  structures  in  the  catalyst  layer  and  multiplied  by  the 
term  (1  —  s)  to  account  for  the  occlusion  of  reaction  sites 
due  to  the  accumulation  of  liquid  water  within  the  cell  [17]. 


The  effectiveness  term,  given  by  Eq.  (25),  is  a  measure 
of  how  readily  reactants  diffuse  through  the  spherical  ag¬ 
glomerates.  An  effectiveness  of  1.0  indicates  that  reactants 
diffusing  through  the  agglomerates  encounter  no  resistance. 
An  effectiveness  less  than  1.0  indicates  that  the  agglomer¬ 
ate  offers  resistance  to  reactant  diffusion  thereby  limiting 
the  reaction  rate.  The  rate  of  reaction  is  controlled  by  the 
reactant  concentration  in  the  polymer  at  the  interface  with 
the  reactant  gas  as  given  by  Henry’s  law,  Eq.  (27),  and 
by  the  local  activation  overpotential  given  by  (0e,d  —  </>i). 
The  electrical  potential  is  assumed  to  be  constant  over  each 
electrode  (catalyst  and  GDL).  As  indicated  by  Eq.  (24),  the 
electrical  potential  is  set  to  zero  on  the  anode  side  and  to 
the  negative  of  the  total  overvoltage  on  the  cathode  side. 

Conservation  of  thermal  energy  is  expressed  by  Eq.  (8) 
and  reflects  heat  transfer  by  conduction  and  convection  as 
well  as  source  terms  for  Ohmic  heating  due  to  ionic  resis¬ 
tance  as  given  by  Eq.  (17),  reversible  heat  as  defined  by 
Eq.  (18),  heat  produced  via  activation  losses  as  given  by 
Eq.  (19),  and  the  latent  heat  associated  with  the  phase  change 
of  water  as  given  by  Eq.  (20).  In  this  model,  the  collector 
plates  are  not  included  in  the  solution  domain.  In  an  actual 
fuel  cell,  the  shoulder  of  the  plates  would  be  in  contact  with 
the  MEA  and  provide  a  low-resistance  pathway  for  heat.  The 
source  term,  5pi,  which  is  given  by  Eq.  (21)  approximates 
the  heat  transfer  through  the  collector  plates  had  they  been 
part  of  the  solution  domain. 

2.2.  Mass  transfer  between  dissolved,  liquid, 
and  vapor  phases 

The  model  presented  here  assumes  that  the  three  phases 
in  which  water  can  exist — liquid,  vapor,  and  dissolved — are 
in  equilibrium.  With  this  assumption,  the  model  does  not 
address  the  actual  rate  of  transport  between  the  three  phases 
or  the  actual  path  the  water  follows  in  moving  among  the 
phases.  Instead,  the  phases  are  simply  assumed  to  be  con¬ 
nected  by  at  least  one  transport  path  with  the  mass  transfer 
coefficients  along  this  path  sufficiently  large  to  ensure  equi¬ 
librium.  The  water  in  the  dissolved  phase  is  assumed  to 
be  exchanged  with  the  vapor  phase.  For  convenience  and 
numerical  stability,  it  is  assumed  that  water  leaving  the 
dissolved  phase  passes  through  the  vapor  phase  and  goes 
directly  into  the  liquid  phase.  Water  entering  the  dissolved 
phase  comes  directly  from  the  vapor  phase.  Water  is  ex¬ 
changed  freely  between  the  liquid  and  vapor  phases,  thus 
maintaining  all  three  phases  in  equilibrium. 

Mass  transfer  among  the  phases  and  equilibrium  con¬ 
straints  are  implemented  by  the  source  terms  of  Eqs.  (14) 
and  (15)  and  switch  functions  given  by  Eqs.  (37)  and  (38). 
The  equilibrium  between  liquid  and  vapor  phases  is  main¬ 
tained  by  Eq.  (14)  which  yields  a  large  positive  value  for 
the  term  S\y  (corresponding  to  rapid  evaporation)  if  liquid 
is  present  and  the  adjacent  vapor  exists  at  a  pressure  below 
saturation.  Conversely,  if  the  pressure  of  the  vapor  phase  ex¬ 
ceeds  the  saturation  pressure,  the  term  S\y  exhibits  a  large 
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Table  4 


Boundary  conditions 


Governing  equation 

Anode  inlet 

Anode  exit 

Cathode  inlet 

Cathode  exit 

Comments 

Momentum 

uy,  avg  —  1547  mm/s 

P  —  P  exit 

uy, avg  =  3747  mm/s 

P  —  P  exit 

— 

Oxygen  transport 

V  •  wo2  =  0 

V  •  wo2  =  0 

wo 2  =  0.21 

V  •  wo2  —  0 

Initial  concentration  in  the  anode 

is  set  to  zero 

Water  vapor  transport 

mjwv  —  0.62 

V  •  uj\vv  —  0 

MJWV  =0.10 

V  •  wj\vv  —  0 

— 

Dissolved  water  transport 

< 

II 

o 

< 

II 

O 

< 

cs 

II 

o 

< 

II 

O 

Initial  concentration  based  on 
equilibrium  with  the  water  vapor 

Liquid  water  transport 

s  =0 

V  •  J  =  0 

s  =0 

V  •  J  =  0 

— 

Membrane  potential 

<1 

II 

o 

V  •  0i  =  0 

<1 

II 

o 

V  •  0i  =  0 

— 

Energy 

353  K 

— 

353  K 

— 

Constant  temperature  BCs  along 
the  gas  channels 

negative  value  (corresponding  to  rapid  condensation).  The 
dissolved  phase  is  maintained  in  equilibrium  by  Eq.  (15) 
which  transfers  water  from  the  dissolved  phase  to  the  va¬ 
por  phase  in  proportion  to  the  difference  between  the  vapor 
density,  p^v  and  the  density  of  the  vapor  in  equilibrium 
with  the  dissolved  phase,  p^v-  Water  evaporating  from  the 
dissolved  phase  passes  directly  through  the  vapor  phase  and 
into  the  liquid  as  reflected  by  the  term  Swd(1  — YwdVpwl 
in  Eq.  (5).  Water  enters  the  dissolved  phase  from  the  vapor 
phase  as  reflected  by  the  term  /wd^wd  in  the  conservation 
equations  for  overall  mass,  Eq.  (1)  and  water  vapor,  Eq.  (4). 
The  variable  ywD  is  a  switch  term,  defined  by  Eq.  (38) 
that  has  a  value  of  1  when  water  is  leaving  the  dissolved 
phase  and  zero  otherwise.  Finally,  water  entering  or  leav¬ 
ing  the  dissolved  phase  is  reflected  by  the  term  Swd/^Av 
in  the  equation  for  dissolved  water  transport,  Eq.  (6). 2  The 
mass  transfer  coefficients  that  permit  exchange  between  the 
phases  and  the  parameters  for  the  gradual  switch  function, 
yLV,  were  chosen  so  that  the  three  phases  remain  near  equi¬ 
librium  with  one  another.  Numerical  tests  were  conducted 
to  verify  that  near-equilibrium  conditions  were  maintained 
among  the  three  phases  and  that  the  vapor  pressure  at  which 
evaporation/condensation  occurred  was  within  ±2.5%  of  the 
saturation  vapor  pressure  at  the  local  temperature. 

2.3.  Boundary  conditions 

Table  4  contains  the  boundary  conditions  and  the  starting 
solution  used  in  the  model.  Since  the  model  is  solved  by  an 
iterative  solution  technique,  the  choice  of  a  starting  solution 
can  affect  convergence  and  solution  time.  The  starting  so¬ 
lution  for  the  species  was  set  equal  to  their  respective  inlet 
boundary  values. 


2  The  term  Swd  has  dimensions  of  mass  of  water  produced  per  unit  time 
per  unit  volume,  g/mm3  s  which  is  consistent  with  gas  species  transport 
Eqs.  (1),  (3),  and  (4).  The  liquid  water  transport  equation  is  expressed  in 
terms  of  fraction  of  the  local  volume  filled  with  water  so  that  Swd  must 
be  divided  by  the  liquid  water  density  to  yield  units  of  mm3 -water/mm3  s 

or  simply  s-1.  The  dissolved  water  transport  equation  is  expressed  in 
terms  of  concentration  so  that  Swd  must  be  divided  by  the  molar  mass 
of  water  to  yield  units  of  mol- water/mm3  s. 


2.4.  Physical  property  evaluation 

In  addition  to  governing  equations  and  boundary  condi¬ 
tions,  the  numerical  model  incorporates  a  number  of  pa¬ 
rameters  some  of  which  are  fundamental  physical  properties 
and  others  which  may  be  chosen  by  the  fuel  cell  designer. 
Table  5  indicates  values  for  properties  such  as  the  viscosity 
of  the  reactant  gases,  the  exchange  current  density  on  bright 
platinum,  etc.  that  are  not  design  variables.  In  general,  the 
values  in  Table  5  were  determined  from  the  literature.  Values 
for  the  exchange  current  density,  transfer  coefficient,  num¬ 
ber  of  electrons  in  the  rate  limiting  step,  solubility  of  oxy¬ 
gen  in  Nation,  and  the  diffusion  coefficient  for  oxygen  in 
Nation  were  estimated  from  data  presented  by  Zhang  et  al. 
[19].  The  estimation  of  these  properties  was  required  given 
that  that  data  in  [19]  is  presented  over  a  temperature  range 
of  303-343  K,  and  the  base  case  model  was  run  at  a  temper¬ 
ature  of  353  K.  In  addition,  the  value  used  for  the  diffusion 
coefficient  of  oxygen  through  Nation  is  about  one  half  of 
that  reported  in  [19].  This  was  done  because  the  diffusion 
coefficient  in  recast  Nation  is  smaller  than  in  Nation  mem¬ 
branes  [5].  Additional  properties  for  Nation  were  taken  from 
Zawodzinski  et  al.  [23]. 

Table  6  indicates  values  for  parameters  that  can  be  spec¬ 
ified  by  the  fuel  cell  designer.  Values  in  Table  6  were  gen¬ 
erally  determined  by  direct  measurements  or  specifications 
for  the  base  case  fuel  cell  that  was  used  for  validation  of  the 
numerical  model.  The  base  case  fuel  cell  used  for  validation 
consisted  of  a  5  cm2  fuel  cell  assembly,  carbon  flow  field 
plates,  ELAT®  GDLs,  and  a  catalyzed  Nation  112  mem¬ 
brane.  Geometric  properties  of  the  fuel  cell  assembly  were 
measured  directly.  In  the  assembled  cell,  GDL  thickness  is 
non-uniform  as  it  is  highly  compressed  under  the  flowfield 
shoulders,  but  not  under  the  gas  channels.  The  thickness  of 
the  GDL  was  calculated  as  the  average  of  the  uncompressed 
GDL  thickness  and  the  compressed  thickness,  which  was  set 
by  the  gasket  material.  The  porosity  of  the  uncompressed 
GDL  was  assumed  to  be  0.6.  This  value  was  then  adjusted 
to  account  for  a  decrease  in  void  volume  due  to  GDL  com¬ 
pression.  The  tortuosity  factor  that  modifies  the  gas  species 
diffusion  coefficients  in  the  GDL  was  estimated  from  the 
work  of  Springer  et  al.  [22]. 
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Table  5 


Physical  properties 


Property 

Value 

Source 

Faraday’s  constant,  F 

96487  C/mol 

[8] 

Permeability  of  gas  diffusion  layer,  /cgdl 

1.8  x  10-5  mm-2 

[8] 

Cathode  gas  viscosity,  /xa ir 

1.0  x  10-5  Pas 

Calc.a 

Anode  gas  viscosity,  /x Ho 

2.0  x  10-5  Pas 

Calc.a 

Liquid  water  viscosity  at  80  °C,  |jiwl 

4.0  x  10_3g/(cms) 

[26] 

Diffusivity  of  oxygen  in  the  polymer, 

2.0  x  10-4mm/s 

Estb 

Diffusivity  of  hydrogen  in  the  polymer,  Dm  H2 

7.9  x  10_4mm/s 

[7] 

Reference  anode  exchange  current  density,  z'o,h2 

3.0  x  10-5  A/mm2 

[27] 

Reference  cathode  exchange  current  density,  io,o2 

4.1  x  10-9A/mm2 

Est.b 

Anodic  transfer  coefficient,  aa 

0.50 

[27] 

Cathodic  transfer  coefficient,  ac 

0.55 

Est.b 

Oxygen  reference  concentration,  co2,ref 

1.18  x  10-9  mol/mm3 

Est.b 

Hydrogen  reference  concentration,  cn2,ref 

2.66  x  10-8  mol/mm3 

Cale.c 

Solubility  coefficient  for  the  cathode,  h^c 

0.19 

Est.b 

Solubility  coefficient  for  the  anode,  /zd,a 

0.64 

Cale.c 

Entropy  of  reaction — anode,  sj? 

42.5  J/(mol  K) 

Calc.d 

Entropy  of  reaction — cathode,  s9 

1 ,  C 

126.8  J/(molK)  (liq  wtr),  64.8  J/mol-K  (wtr  vpr) 

Calc.d 

a  Calculated  from  inlet  conditions. 

b  Estimated  from  data  in  [19]  at  a  temperature  of  343  K  and  fully  humidified  O2  pressure  of  1  atm. 
c  Calculated  from  data  in  [7]  at  a  temperature  of  353  K  and  H2  pressure  of  1  atm. 
d  Based  on  anode  and  cathode  half  reactions. 


The  catalyst  layer  for  the  base  case  cell  was  prepared 
using  a  catalyst  ink  composed  of  carbon  supported  cata¬ 
lyst  (20  wt.%  Pt  on  Vulcan  XC-72R)  dispersed  in  a  5  wt.% 
Nation®  1100  solution.  The  fabrication  technique  was  sim¬ 
ilar  to  the  decal  method  given  in  Wilson  et  al.  [20]  with 
the  exception  that  the  ink  was  applied  in  the  protonated,  not 
TBA+  form.  The  mass  and  composition  of  the  catalyst  layer 
were  recorded  and  used  to  determine  the  volume  fraction  of 
the  various  constituents  and  the  overall  catalyst  layer  poros¬ 
ity.  The  thickness  of  the  catalyst  layer  was  determined  from 
SEM  images  of  the  MEA  cross-section  such  as  the  one  pre¬ 
sented  in  Fig.  2. 

The  catalyst  layer  porosity  was  determined  by  evaluating 
the  solid  volume  of  the  catalyst  layer  based  on  its  mass 
and  composition  and  the  total  volume  of  the  catalyst  layer 
based  on  its  thickness.  The  mass  of  the  catalyst  layer  was 


Fig.  2.  SEM  image  of  the  test  MEA.  Magnification  is  5000x. 


determined  by  weighing  the  catalyst  layer  after  it  had  been 
painted  onto  a  Teflon  decal  and  dried  to  remove  all  of  the 
solvent.  The  solid  volume  of  the  catalyst  layer,  Vs,  was  then 
be  determined  by 


^cat-Xp  mca  tvCTpt  mcatv  c(l 

Pp  PPt  PC 


(40) 


where  mcat  is  the  mass  of  the  catalyst  layer,  vp  the  mass 
fraction  of  polymer  in  the  catalyst,  xq  the  mass  fraction  of 
carbon  supported  catalyst,  vpt  the  mass  fraction  of  platinum 
in  the  carbon  supported  catalyst,  pv  the  density  of  the  poly¬ 
mer,  pc  the  density  of  the  carbon  support,  yOpt  the  density 
of  platinum.  The  total  volume  of  the  catalyst  layer  includ¬ 
ing  the  pores,  Vtot,  can  be  estimated  from  the  SEM  images. 
From  these  two  volumes  the  porosity,  e^d,  which  is  the 
volume  fraction  of  pores  in  the  catalyst  layer  relative  to  the 
total  layer  volume,  can  be  calculated  and  is  expressed  as 


Ccat 

£void 


Vtot  ~  Es 

Vtot 


(41) 


This  approach  for  determining  porosity  assumes  that  the 
Nation  and  carbon  phases  do  not  coexist  (i.e.  the  Nation 
does  not  fill  the  micropores  in  the  activated  carbon).  This 
assumption  is  consistent  with  the  work  of  Uchida  et  al.  [24] 
who  concluded  that  the  Nation  only  exists  in  the  larger  pores 
such  as  those  between  particles  or  agglomerates. 

Cyclic  voltammetry  (C V)  was  used  to  determine  the  active 
catalyst  area,  Aact,  within  the  catalyst  layer.  The  CV  tests 
were  conducted  using  an  approach  similar  to  that  described 
by  O’Hayre  et  al.  [25].  With  this  approach,  the  electrical 
potential  applied  to  the  catalyst  layer  is  varied  from  —0.1  to 
1.0  V  and  back  in  a  triangle  waveform  while  the  current  is 
recorded.  The  current  is  integrated  over  time  to  determine 
the  charge  transferred  during  the  adsorption  and  desorption 
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Table  6 


Design  and  operating  parameters  for  validation  study 


Property 

Value 

Source 

Gas  channel  width,  Wgc 

30.0  cm 

Measured 

Gas  channel  length,  Lgc 

1.0  mm 

Measured 

Gas  channel  height,  Hgc 

1.0  mm 

Measured 

Collector  thickness,  tco\ 

1.0  mm 

Measured 

Anode  GDL  thickness,  tagdi 

0.290  mm 

Est.a 

Cathode  GDL  thickness,  fcgdi 

0.254 

Est.a 

Gas  diffusion  layer  void  fraction,  £^jd 

0.375 

Est.a 

Catalyst  layer  thickness,  tcat 

0.0165  mm 

Measured 

Pt  and  carbon  volume  fraction  in  the 

0.45 

Calc.b 

catalyst  layer,  sfal 

Catalyst  layer  void  fraction,  £^atid 

0.31 

Calc.b 

Polymer  volume  fraction  in  the  catalyst 

0.24 

Calc.b 

layer,  £pat 

Membrane  thickness,  tm 

0.0508  mm 

Measured 

Cell  temperature,  Tco\\ 

353  K 

Measured 

Inlet  pressure,  Pm 

310kPa 

Measured 

Air  inlet  relative  humidity,  Rhc 

100% 

Measured 

Fuel  inlet  relative  humidity,  Rha 

100% 

Measured 

Theoretical  open  circuit  voltage,  E{ h 

1.19 

Calc.c 

Open  circuit  voltage,  Vqc 

Varies 

Measured 

Specific  reaction  area  of  the  catalyst 

6990  mm-1 

Measured 

layer,  Av  i 

Mean  agglomerate  size,  Laggi 

400.0  nm 

Data  fitd 

Tortuosity  of  the  GDL,  Tgdl 

3.5 

Data  fitd 

Evap./cond.  mass  transfer  coefficient,  f 

2  g/mm3  s 

Equil.e 

Dissolved/vapor  mass  transfer 

5000s-1 

Equil.f 

coefficient,  hm 

a  Estimated  from  the  uncompressed  thickness,  void  fraction,  and  degree 
of  compression. 

b  Calculated  from  the  catalyst  layer  composition  and  mass. 

c  Calculated  from  Nernst  equation  for  base  case  reactant  temperature, 
pressure,  and  composition. 

d  Property  adjusted  (within  the  range  reported  in  the  literature)  to  fit 
the  data. 

e  Chosen  large  enough  to  maintain  equilibrium  between  liquid  and 
vapor  phases. 

f  Chosen  large  enough  to  maintain  equilibrium  between  dissolved  and 
vapor  phases. 


of  a  monolayer  of  hydrogen  on  the  active  catalyst  surface. 
The  active  area  of  the  catalyst  can  be  then  be  calculated  from 


A 


act  — 


Qa 

0pt 


(42) 


diffusive  flux  of  the  gas  species  through  the  membrane, 
which  is  assumed  impermeable  to  gases,  the  diffusivities  of 
all  gas  species  are  set  to  zero  within  the  membrane  region 
of  the  model.  The  domain  is  divided  into  64  x  121  elements. 
Mapped  meshing  is  used  to  maintain  a  sufficient  mesh 
density  throughout  the  model  domain.  In  regions  of  large 
gradients,  such  as  at  interfaces  and  in  the  catalyst  layers, 
mesh  size  is  decreased,  while  in  the  gas  channels  and  other 
regions  of  relatively  small  gradients,  it  is  much  coarser.  A 
sensitivity  analysis  was  conducted  by  doubling  the  number 
of  elements  in  the  mesh.  The  solution  changed  on  average 
by  less  than  2.0%  and  so  was  assumed  to  be  mesh  indepen¬ 
dent.  The  equations  are  solved  with  the  commercial  CFD 
solver,  CFDesign®. 


3.  Results  and  discussion 

The  numerical  model  is  validated  by  comparing  model 
results  to  experimental  data.  The  base  case  cell  was  run  at 
a  cell  temperature  of  80  °C  with  the  inlet  reactant  gases  on 
both  the  anode  and  cathode  sides  maintained  at  a  tempera¬ 
ture  of  80  °C,  a  pressure  of  30psig,  and  a  relative  humidity 
of  100%.  The  mass  flow  rate,  at  both  the  anode  and  cath¬ 
ode,  corresponded  to  a  stoichiometric  ratio  of  6  at  a  current 
density  of  1  A/cm2.  A  comparison  of  model  results  with  ex¬ 
perimental  data  at  30  psig  (base  case)  and  20  psig  is  shown 
in  Fig.  3a  and  b,  respectively. 


where  Qa  is  the  amount  of  charge  (jxC)  transferred  during 
adsorption  (or  desorption)  and  <2pt  the  charge  transferred 
during  the  adsorption  of  a  monolayer  of  hydrogen  on  atom¬ 
ically  smooth  platinum  (210  puC/cm2). 

2.5.  Numerical  methods 

The  outer  surfaces  of  the  gas  channels  shown  in  Fig.  1 
are  bounded  by  the  collector  plates  (not  shown)  and  are 
impermeable  to  gases.  As  part  of  the  single  domain  for¬ 
mulation,  each  governing  equation  is  solved  throughout  the 
entire  domain,  even  if  the  equation  is  not  physically  valid 
in  every  region.  This  is  accomplished  using  a  variety  of 
numerical  techniques  [8,21].  For  instance,  to  eliminate  the 


Fig.  3.  (a)  Model  comparison  with  test  data  at  30 psig  (base  case),  (b) 
Model  comparison  with  test  data  at  20  psig. 
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Fig.  4.  Effect  of  liquid  water  buildup  on  cell  performance. 


Some  of  the  modeling  parameters  were  estimated  by  ad¬ 
justing  their  values  until  good  agreement  was  achieved  be¬ 
tween  the  test  results  and  the  model  output.  This  was  done 
only  once,  for  the  30psig  case.  The  results  shown  in  Fig. 
3b  for  20  psig  indicate  that,  without  further  adjustment,  the 
model  is  able  to  accurately  predict  performance  when  the 
cell  is  running  under  conditions  that  differ  from  the  base 
case. 

Of  particular  interest  in  Figs.  3a  and  b  is  the  loss  in  per¬ 
formance  due  to  mass  transfer  limitation  at  high  current  den¬ 
sity.  The  increasingly  steep  drop  in  performance  as  the  cell 
approaches  its  limiting  current  density  is  attributable  to  the 
saturation  of  the  cathode  GDL  and  catalyst  layer  with  wa¬ 
ter  and  the  corresponding  restriction  in  reactant  transport. 
Previous  efforts  [3,12]  to  explain  limiting  current  behavior 
based  solely  on  diffusive  resistance  in  the  agglomerates  led 
to  a  more  abrupt  drop  in  performance.  By  including  the  ef¬ 
fects  of  liquid  water  transport,  the  current  model  is  able  to 
more  closely  simulate  performance  in  the  region  where  mass 
transport  effects  begin  to  dominate. 

The  importance  of  liquid  water  transport  to  the  accurate 
modeling  of  fuel  cell  performance  is  further  illustrated  in 
Figs.  4  and  5.  Performance  curves  with  and  without  liquid 
water  transport  are  shown  for  the  base  case  conditions  in 
Fig.  4.  The  curve  labeled  “NoSat”  does  not  include  the  ef¬ 
fect  of  liquid  water  accumulation  on  gaseous  reactant  trans¬ 
port.  The  predicted  performance  when  liquid  water  effects 
are  neglected  is  much  better  than  for  the  case  labeled  “Sat”, 


Fig.  5.  Influence  of  liquid  water  on  current  density  variation  along  the 
flow  direction. 


which  does  include  the  effect  of  liquid  water  accumulation 
on  reactant  transport.  Comparison  of  the  two  curves  demon¬ 
strates  that  the  effects  of  liquid  water  accumulation  become 
apparent  even  at  relatively  low  values  of  current  density. 
Furthermore,  when  liquid  water  effects  are  not  included  in 
the  model,  the  cell  voltage  does  not  exhibit  an  increasingly 
steep  drop  as  the  cell  approaches  its  limiting  current  density. 
This  drop  off  in  performance  is  clearly  demonstrated  by  ex¬ 
perimental  data,  but  cannot  be  accurately  modeled  without 
the  incorporation  of  liquid  water  transport. 

The  variation  of  current  density  along  the  channel  for 
the  same  two  cases  at  a  cell  voltage  of  0.4  V  is  shown  in 
Fig.  5.  In  the  case  where  liquid  water  does  not  impede  re¬ 
actant  transport  (“NoSat”)  the  current  density  produced  in 
the  cell  is  greater  in  magnitude  down  the  channel  relative  to 
the  more  physically  realistic  case  (“Sat”).  In  the  absence  of 
liquid  water,  the  drop  in  current  density  as  the  reactants  are 
depleted  down  the  channel  is  more  pronounced.  When  satu¬ 
ration  is  considered,  a  decrease  in  current  density  in  the  flow 
direction  reduces  the  water  production  rate  and  leads  to  a 
lower  degree  of  saturation  of  the  cathode  GDL  and  catalyst 
layer  at  the  exit  to  the  cell.  With  less  accumulated  water,  the 
GDL  and  catalyst  layer  are  more  open  to  reactant  flow,  thus 
counteracting  the  effect  of  reactant  depletion  and  leading  to 
more  uniform  current  density  along  the  gas  channel. 

Fig.  6  shows  a  contour  plot  of  liquid  water  saturation  for 
the  base  case  conditions  in  the  cathode  GDL  and  catalyst 
layer  at  an  average  current  density  of  1 . 18  A/cm2.  The  model 
predicts  that  the  level  of  saturation  decreases  along  the  flow 
direction.  The  decrease  in  saturation  is  due  to  the  reduced 
reaction  rate  (associated  with  the  drop  in  reactant  concen¬ 
tration  along  the  flow  direction)  and  a  corresponding  drop 
in  the  amount  of  water  produced  as  well  as  the  amount  of 
water  dragged  across  the  membrane  from  the  anode. 

Results  from  the  model  also  demonstrate  the  signifi¬ 
cance  of  water  transport  from  the  anode  to  the  cathode  by 
electro-osmotic  drag.  This  water  transport  reduces  the  per¬ 
formance  of  the  cell  both  by  dehydrating  the  anode,  thereby 
decreasing  its  conductivity,  and  by  increasing  the  accumu¬ 
lation  of  liquid  water  at  the  cathode,  which  decreases  its 
permeability  to  reactant  gas  flow. 

Fig.  7  shows  profiles  for  polymer  water  content  in  the 
catalyst  layers  and  membrane  as  a  function  of  the  local 
current  density  and  position  through  the  thickness  of  the 
MEA  at  a  point  halfway  down  the  channel  for  the  base  case 
conditions.  The  influence  of  electro-osmotic  drag  is  readily 
apparent  from  these  results.  At  low  current  density,  there 
is  very  little  change  in  water  content  across  the  MEA.  This 
is  due  to  a  relatively  low  amount  of  drag  and  to  the  fact 
that  the  vapor  activity  at  the  anode  and  cathode  is  nearly 
identical.  As  current  density  is  increased,  the  water  con¬ 
tent  profile  becomes  steeper  as  the  anode  dehydrates  and 
the  cathode  water  content  increases.  The  results  shown  in 
Fig.  7  also  indicate  that  as  current  density  increases,  the  to¬ 
tal  amount  of  water  contained  in  the  MEA  decreases.  This 
occurs  because  the  vapor  activity  of  the  anode  stream  has 
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Gas  channel  _  FLOW  CHANNEL  _ fe.  Gas  channel 

inlet  ^  exit 


Fig.  6.  Saturation  profiles  in  the  cathode. 


Fig.  7.  Variation  of  polymer  water  content  through  the  MEA. 


dropped  due  to  water  removal  upstream  leaving  less  water 
available  to  hydrate  the  anode. 

The  effect  of  the  water  content  gradient  as  well  as  the  sig¬ 
nificance  of  the  electro-osmotic  drag  are  further  illustrated 
by  the  data  in  Table  7.  Considering  electro-osmotic  drag 
alone,  the  amount  of  water  moved  from  the  anode  to  the  cath¬ 
ode  for  fully  hydrated  Nation®  would  be  between  2.0  and 
2.9  water  molecules  per  proton  at  30  °C  [23].  Due  to  opera¬ 
tion  with  a  partially  hydrated  membrane  and  catalyst  layers 
as  well  as  back  diffusion  of  water  from  the  cathode  to  the 

Table  7 


Water  transport  through  the  MEA 


Average  current 
density  (A/cm2) 

Net  water  transport 
per  proton 
(H20/H+) 

Fraction  of  water  accumulation 
at  the  cathode  due  to  transport 
across  the  MEA 

0.05 

0.16 

0.24 

0.52 

0.34 

0.41 

1.18 

0.30 

0.37 

anode,  the  net  amount  of  water  moved  per  proton  predicted 
by  the  model  is  between  0.16  and  0.34,  which  is  consistent 
with  the  modeling  results  presented  by  others  [6,7].  It  is  in¬ 
teresting  to  note  that  the  water  transported  across  the  MEA, 
from  the  anode  to  the  cathode,  makes  up  between  20  and 
40%  of  the  total  amount  of  water  accumulation  at  the  cath¬ 
ode  (transport  and  electrochemical  production).  This  indi¬ 
cates  that  in  addition  to  limiting  cell  performance  by  way  of 
anode  dehydration,  water  transport  by  electro-osmotic  drag 
is  responsible  for  a  significant  fraction  of  the  liquid  water 
buildup  at  the  cathode  and  the  resulting  reactant  transport 
limitations. 


4.  Conclusions 

A  comprehensive,  steady  state,  two-dimensional  model 
including  liquid  water  transport  has  been  developed  and 
validated  with  experimental  data.  Results  from  the  model 
show  that,  in  order  to  accurately  simulate  fuel  cell  op¬ 
eration,  liquid  water  transport  within  the  cell  must  be 
considered  as  it  results  in  a  loss  of  performance  even  at 
relatively  low  current  density.  In  addition,  water  transport 
through  the  polymer  portion  of  the  catalyst  layer  and  the 
membrane  plays  an  important  role  with  respect  to  both 
Ohmic  losses  and  reactant  transport  restrictions  at  the  cath¬ 
ode.  The  model  predicts  a  net  amount  of  water  transport 
across  the  membrane  of  between  0.16  and  0.34  mole  of 
water  per  mole  of  protons  transported  from  the  anode  to  the 
cathode.  This  accounts  for  20-40%  of  the  total  amount  of 
water  accumulation  at  the  cathode,  which  is  a  combination 
of  water  produced  electrochemically  and  water  transported 
by  electro-osmotic  drag.  When  operating  at  high  inlet 
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humidity,  any  water  accumulating  at  the  cathode  remains  in 
the  liquid  phase  and  fills  the  porous  regions  of  the  catalyst 
layer  and  the  GDL.  The  fraction  of  pores  filled  with  liquid 
water  is  highest  within  the  catalyst  layer  near  the  inlet  to 
the  cell. 
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